function [beta_vary,beta_b_vary] = beta_variation(divergent_fraction,divergent_fraction_b,...
        beta,beta_b,covariates,round_digits,filename,rounding_calculations)

% Function which assesses the variation of the beta values obtained by the
% end of the specified number of iterations (if all the non-NaNs are the
% same value, of course, that would mean that there is good robustness in
% the process)

% Measures of the variation of beta is comparison to the mean
beta_vary = nanstd(beta) ./ nanmean(beta);
beta_b_vary = nanstd(beta_b) ./ nanmean(beta_b);
beta_max = max(max(beta_vary),max(beta_b_vary));
beta_min = min(min(beta_vary),min(beta_b_vary));
beta_range = beta_max - beta_min;
% Converting to log scale
beta_vary_log = log10(abs(beta_vary));
beta_b_vary_log = log10(abs(beta_b_vary));
beta_log_max = max(max(beta_vary_log),max(beta_b_vary_log));
beta_log_min = min(min(beta_vary_log),min(beta_b_vary_log));
beta_log_range = beta_log_max - beta_log_min;


% Discover the names of the covariates chosen
% NOTE: Explicitly done for specific data files
if (strcmp(filename,'myopia_edited.csv'))
    all_covariate_names = ["age", ...
        "spheq", ...
        "al", ...
        "acd", ...
        "it", ...
        "vcd", ...
        "sporthr", ...
        "readhr", ...
        "comphr", ...
        "studyhr", ...
        "tvhr", ...
        "diopterhr", ...
        "mommy", ...
        "dadmy", ...
        "studyyear"];
elseif (strcmp(filename,'polypharm_edited.csv'))
    all_covariate_names = ["cormorbid", ...
        "anyprim", ...
        "numprim", ...
        "age", ...
        "intptmhv3", ...
        "year", ...
        "urban", ...
        "male", ...
        "white", ...
        "hispanic"];
end


for k=1: numel(covariates)
    covariate_names(k) = all_covariate_names(covariates(k));
end

% Plot the variations on standard scale 
figure(500)
hold on
plot(beta_vary,'ko','MarkerSize',7)
plot(beta_b_vary,'kx','MarkerSize',8)
xlim([0 numel(covariates)+1])
%xlabel('covariates $ j $','Interpreter','Latex')
ylabel('$ \sigma(\beta_j) / \bar{\beta}_j $', 'Interpreter','Latex')
xtickangle(35)
set(gca,'xtick',[1:numel(covariates)],'xticklabel',covariate_names)
if(strcmp(rounding_calculations,'yes'))
    title(['Round digits = ',num2str(round_digits)])
end
text(numel(covariates)-2,beta_min+0.98*beta_range,...
            ['SAMLE:  ',num2str(divergent_fraction)],'fontsize',14)
text(numel(covariates)-2,beta_min+0.92*beta_range,...
            ['non-SAMLE:  ',num2str(divergent_fraction_b)],'fontsize',14)
ax = gca;
ax.XLabel.FontSize = 30; % This one doesn't seem to work, for some reason (seems to want ticklabels to be same fontsize as this)
ax.YLabel.FontSize = 18;
ax.XAxis.FontSize = 12; %XTickLabel fontsize 
                       %(do this rather than mess with XTickLabel, because
                       %these are nonnumeric)
ax.YAxis.FontSize = 18;
ax.Title.FontSize = 18;

% Plot the variations on a log scale to be able to visualize the higher
% variations while also seeing variations of small size
figure(501)
hold on
plot(beta_vary_log,'ko','MarkerSize',7)
plot(beta_b_vary_log,'kx','MarkerSize',8)
xlim([0 numel(covariates)+1])
%xlabel('covariates $ j $','Interpreter','Latex')
ylabel('$ \log_{10} \left| \sigma(\beta_j) / \bar{\beta}_j \right|$', 'Interpreter','Latex')
xtickangle(35)
set(gca,'xtick',[1:numel(covariates)],'xticklabel',covariate_names)
if(strcmp(rounding_calculations,'yes'))
    title(['Round digits = ',num2str(round_digits)])
end
text(numel(covariates)-2,beta_log_min+0.98*beta_log_range,...
            ['SAMLE:  ',num2str(divergent_fraction)],'fontsize',14)
text(numel(covariates)-2,beta_log_min+0.92*beta_log_range,...
            ['non-SAMLE:  ',num2str(divergent_fraction_b)],'fontsize',14)
ax = gca;
ax.XLabel.FontSize = 30; % This one doesn't seem to work, for some reason (seems to want ticklabels to be same fontsize as this)
ax.YLabel.FontSize = 18;
ax.XAxis.FontSize = 12; %XTickLabel fontsize 
                       %(do this rather than mess with XTickLabel, because
                       %these are nonnumeric)
ax.YAxis.FontSize = 18;
ax.Title.FontSize = 18;

